Tomography, the computation (estimation) of densities (represented as voxels) in a region of n-dimensional space based on m-dimensional projections of that region (represented as pixels; usually 0<m<n), falls into two major categories: linear filtered back-projection (FBP) and fast Fourier transform (FFT) methods and modern non-linear iterative reconstruction (IR), such as the Algebraic Reconstruction Technique (ART), the Simultaneous Iterative Algebraic Reconstruction Technique (SIRT and SART), and the Model Based Iterative Reconstruction (MBIR). Image reconstruction from modern tomography settings, however, is typically complex due to the use of non-linear estimation. Tomographic settings often require 106-1010 image elements (i.e., voxels) to be computed when one projection contains anywhere from 2,000 pixels (one-dimensional) to 4,000×4,000 pixels (two-dimensional).
Major deficiencies associated with FBP include the need for a large number of projections to achieve limited quantitative accuracy. The number of projections is typically counted in the hundreds or thousands but the projections are not used as efficiently as they could be.
The benefits of IR techniques include their ability to reduce reconstruction errors when following FBP. An overview of related processing techniques and their particular benefits, especially linear and robust techniques, is presented in Sunnegardh [15]. Additionally, IR can account for constraints, especially the ability to assure density estimates to be non-negative. Finding the optimal object reconstruction requires the minimization of the objective function, which is the sum of all reconstruction errors. Finding the optimal reconstruction requires operations on each voxel. Near the minimum error the objective function can typically be represented by the Hessian matrix, describing the second derivative of the objective function. In tomography, for a three-dimensional volume of 1,000 voxels or more in any one dimension, equivalent to 1,000×1,000×1,000=109 voxels or more, this Hessian may contain 109×109=1018 or more elements. Storing and manipulating such a large number of elements, however, is out of reach for present computers, and for systems in the foreseeable future.
To date, all commercial IR systems performance has been restrained by the properties of this Hessian matrix. Iterations, using quasi-Newton methods, bypass the evaluation of the Hessian matrix. For these methods, the tendency is to identify minimization steps that cope with the largest eigenvalues of the Hessian. For this purpose, an initial FBP is followed up by iterative refinement. Due to the large size of the Hessian matrix for image reconstruction, however, its structure (and that of its inverse) are typically ignored or poorly approximated during the refinement. Furthermore, because of the wide distribution of eigen-values of the Hessian, current optimization techniques tend to show no improvement beyond a number of iterations (typically counted in the tens-to-thousands) and may only cope with few of the large eigenvalues.
Multi-grid variations of these algorithms may help, but ultimately still fail because of the size of the Hessians involved with fine grids. Multi-grid resolution here refers to the use of progressively finer resolution as iterations are performed.
Further, systems uncertainties in tens or hundreds of systems parameters, characterizing, for example, object deformation, beam hardening and scatter in the case of x-ray imaging, or field distortions in the case of magnetic resonance imaging (MRI), may also increase measurement data inconsistencies.